//
//  CFDProblem.h
//  NSESolver
//
//  Created by Patrick Zulian on 3/29/11.
//  Copyright 2011 ICS - USI - Lugano. All rights reserved.
//

#ifndef NSESolver_CFDProblem_H
#define NSESolver_CFDProblem_H

#include "Types.h"
#include "BoundaryHandler.h"
#include "InitialValue.h"
#include "CartesianTopology.h"

namespace nsesolver{
    /**
     * \brief The CFDProblem class ... TODO Doc
     */ 
    class CFDProblem {
    private:
        typedef mpi::CartesianTopology CartesianTopology;
    
        const int imax, jmax, iGlobal, jGlobal;
        const double dx, dy, alpha;
        
        Vector2d southWestCorner;
        Vector2d northEastCorner;
        
        BoundaryHandler bc;
        InitialValue iv;
        
        Matrixd u, v, p;
        
        CartesianTopology cart;
        
        
        static int computeMaxLoc(const CartesianTopology &cart, const int amax, const int dim);
        static int computeGlobalIndex(const CartesianTopology & cart, const int amax, const int dim);
    public:
        
       /**
        * \brief Constructs a CFD Problem arranging the decomposition of the domain according to the cartesian topology
        * \param imax the global max i index [0, imax] is the index range
        * \param jmax the global max j index [0, jamx] is the index range
        * \param southWestCorner the x,y values of the corner
        * \param northEastCorner the x,y values of the corner
        * \param bc is the object the manages all the boundary conditions and in the parallel case also the ghosts nodes communication
        * \param iv initalizes u0 and v0
        * \param cart the mpi based cartesian topology
        */
        CFDProblem(const int imax, const int jmax, const double alpha,
                   const Vector2d & southWestCorner, const Vector2d & northEastCorner, 
                   const BoundaryHandler &bc, const InitialValue &iv, const CartesianTopology &cart);
        
        /**
         * The Destructor
         */
        ~CFDProblem();
        
        
        /// \return the topology of the processeor grid associated with this problem
        inline CartesianTopology & getTopology()
        {
            return cart;
        }
        
         /// \return the topology of the processeor grid associated with this problem
        inline const CartesianTopology & getTopology() const
        {
            return cart;
        }
        
        inline const BoundaryHandler &getBoundaryConditions() const 
        {
            return bc;
        }
        
        inline BoundaryHandler &getBoundaryConditions() 
        {
            return bc;
        }
        
       /**
        * \return the reference to the horizontal velocity
        */
        inline 
        Matrixd & U()
        {
            return u;
        }
        
        /**
         * \return the reference to the horizontal velocity field
         */
        inline 
        const Matrixd & U() const
        {
            return u;
        }
        
        
        /**
         * \return the reference to the vertical velocity field
         */
        inline 
        Matrixd & V()
        {
            return v;
        }
        
        /**
         * \return the reference to the vertical velocity field
         */
        inline 
        const Matrixd & V() const
        {
            return v;
        }
        
        
        /**
         * \return the reference to the pressure field
         */
        inline 
        Matrixd & P()
        {
            return p;
        }
        
        /**
         * \return the reference to the pressure field
         */
        inline 
        const Matrixd & P() const
        {
            return p;
        }
        
        
        /**
         * \return max abs value of the horizontal velocity u
         */
        inline double uAbsMax() const
        {
            double uMax = u.getAbsMax();
            cart.getCommunicator().allReduce(&uMax, 1, MPI_MAX);
            return uMax;
        }
        
        /**
         * \return max abs value of the vertical velocity v
         */
        inline double vAbsMax() const
        {
            double vMax = v.getAbsMax();
            cart.getCommunicator().allReduce(&vMax, 1, MPI_MAX);
            return vMax;
        }
        
        
        /// sets the boundary according to its boundary conditions and domain decomposition for u and v
        inline void setBoundary()
        {
            bc.setBoundary(*this);
        }
        
        /// sets the boundary according to its boundary conditions and domain decomposition for p
        inline void setPBoundary()
        {
            bc.setPBoundary(*this);
        }
        
        
        inline void setFGBoundary(Matrixd &F, Matrixd &G)
        {
           bc.setFGBoundary(F, G, *this);
        }
        
        /// sets the initial values of u and v
        inline void setInitialValues()
        {
            iv.setInitialValues(*this);
        }
        
        /// computes the differentials at point i j
        inline double d2udx2(const int i,const int j) const
        {
            return (u(i+1,j)-2.0*u(i,j)+u(i-1,j))/(dx*dx);
        }
        
        /// computes the differentials at point i j
        inline double d2udy2(const int i,const int j) const
        {
            return (u(i,j+1)-2.0*u(i,j)+u(i,j-1))/(dy*dy);
        }
        
         /// computes the differentials at point i j
        inline double duSquaredx(const int i, const int j) const
        {
            return 
            1.0/(dx*4.0)*(
                          Math::Pow((u(i,j)+u(i+1,j)), (int)2)-
                          Math::Pow((u(i-1,j)+u(i,j)), (int)2)+
                          alpha*(
                                 Math::abs(u(i,j)+u(i+1,j))*
                                 (u(i,j)-u(i+1,j))-
                                 Math::abs(u(i-1,j)+u(i,j))*
                                 (u(i-1,j)-u(i,j))
                                 )
                          );
        }
        
         /// computes the differentials at point i j
        inline double duvdy(const int i, const int j) const
        {
            return
            1.0/(dy*4.0)*(
						  (v(i,j)+v(i+1,j))*
						  (u(i,j)+u(i,j+1)) -
						  (v(i,j-1)+v(i+1,j-1))*
						  (u(i,j-1)+u(i,j))+
                          alpha*(
                                 Math::abs(v(i,j)+v(i+1,j))*
                                 (u(i,j)-u(i,j+1))-
                                 Math::abs(v(i,j-1)+v(i+1,j-1))*
                                 (u(i,j-1)-u(i,j))
                                 )
                          );
            
        }
        
         /// computes the differentials at point i j
        inline double d2vdx2(const int i, const int j) const
        {
            return 1.0/(dx*dx)*(
                                v(i+1,j)-2.0*v(i,j)+v(i-1,j));
        }
        
         /// computes the differentials at point i j
        inline double d2vdy2(const int i, const int j) const
        {
            return 1.0/(dy*dy)*(
                                v(i,j+1)-2.0*v(i,j)+v(i,j-1));
        }
        
         /// computes the differentials at point i j
        inline double dvSquaredy(const int i, const int j) const
        {
            return 1.0/(dy*4.0)*(
                                 Math::Pow(v(i,j)+v(i,j+1), (int)2)-
                                 Math::Pow(v(i,j-1)+v(i,j), (int)2)+
                                 alpha*(
                                        Math::abs(v(i,j)+v(i,j+1))*
                                        (v(i,j)-v(i,j+1))-
                                        Math::abs(v(i,j-1)+v(i,j))*
                                        (v(i,j-1)-v(i,j))
                                        )
                                 );
            
        }
        
         /// computes the differentials at point i j
        inline double duvdx(const int i, const int j) const
        {
            return 1.0/(dx*4.0)*(
                                 (u(i,j)+u(i,j+1))*
                                 (v(i,j)+v(i+1,j))-
                                 (u(i-1,j)+u(i-1,j+1))*
                                 (v(i-1,j)+v(i,j))+
                                 alpha*(
                                        Math::abs(u(i, j)+u(i,j+1))*
                                        (v(i,j)-v(i+1,j))-
                                        Math::abs(u(i-1,j)+u(i-1,j+1))*
                                        (v(i-1,j)-v(i,j))
                                        )
                                 );
        }
        
         /// computes the differentials at point i j
        inline double d2pdx2(const int i, const int j) const
        {
            return 1.0/Math::Pow(dx, (int)2)*(
                                              p(i+1,j)-2.0*p(i,j)+p(i-1,j));
        }
        
         /// computes the differentials at point i j
        inline double d2pdy2(const int i, const int j) const
        {
            return 1.0/Math::Pow(dy, (int)2) *
            ( p(i,j+1)-2.0*p(i,j)+p(i,j-1));
        }
        
         /// computes the differentials at point i j
        inline double dpdx(const int i, const int j) const
        {
            return (p(i+1,j)-p(i,j))/dx;
        }
        
         /// computes the differentials at point i j
        inline double dpdy(const int i, const int j) const
        {
            return (p(i,j+1)-p(i,j))/dy;
        }
        
        
        /// \return the local imax 
        inline 
        int IMax() const
        {
            return imax;
        }
        
        /// \return the local jmax
        inline 
        int JMax() const
        {
            return jmax;
        }
       
        /// \return the global i index for this processor
        inline 
        int IGlobal() const
        {
            return iGlobal;
        }
        
          /// \return the global j index for this processor
        inline 
        int JGlobal() const
        {
            return jGlobal;
        }


        /// \return the delta x
        inline
        double DX() const
        {
            return dx;
        }
        
        /// \return the delta y
        inline 
        double DY() const
        {
            return dy;
        }
        
        
    };
    
}
#endif //NSESolver_CFDProblem_H
